% osl_single_subject_maps.m
%
% Syntax: oil = osl_single_subject_maps(oil);
%
% Function to estimate subject specific COPE and T-statistic maps from 
% group ICA decompositions. OIL stages 1-4 must have been run.
%
% Maps are saved in [oil.source_recon.dirname '/' oil.enveloping.name '/' oil.concat_subs.name '/' oil.ica.name '/' oil.ica_first_level.name].
%
% Maps names are stored in oil.ica_first_level.results.cope_files & oil.ica_first_level.results.tstat_files.
%
% Henry Luckhoo (henry.luckhoo@trinity.ox.ac.uk)
%
% Version 1.4
% 040213

function oil = osl_single_subject_maps(oil)

% Setup stats directory to save the results into
save_dir = [oil.source_recon.dirname '/' oil.enveloping.name '/' oil.concat_subs.name '/' oil.ica.name '/' oil.ica_first_level.name];
if ~isdir(save_dir), mkdir(save_dir);end;

Nics = oil.ica.num_ics;
Nvoxels = size(nii_quickread([oil.source_recon.dirname '/' oil.enveloping.name '/' oil.enveloping.results.source_space_envelopes_results_fnames{oil.concat_subs.sessions_to_do(1)}], oil.enveloping.gridstep),1);
subj_ind=oil.concat_subs.results.subj_ind;
Nsubs=length(subj_ind)-1;

con=[]; for i=1:Nics; con{i}=(1:Nics==i)'; end  % set-up design matrix (identity).


if isfield(oil.ica_first_level, 'spatial_basis_set');
     
    [~,~,scales,~,~] = read_avw(oil.ica_first_level.spatial_basis_set);
    if scales(1) ~= oil.enveloping.gridstep;
        oil.ica_first_level.spatial_basis_set = osl_resample_nii(oil.ica_first_level.spatial_basis_set, [oil.ica_first_level.spatial_basis_set '_ds' num2str(oil.enveloping.gridstep) 'mm'], oil.enveloping.gridstep);
    end
    spat_bas = nii_quickread(oil.ica_first_level.spatial_basis_set,oil.enveloping.gridstep);
    
elseif strcmp(oil.ica.temp_or_spat,'spatial')
    spat_bas = oil.ica.results.sICs';
else
    spat_bas = oil.ica.results.mixing_matrix;
end

if oil.ica_first_level.use_nonorm_data
    concat_data = nii_quickread([oil.source_recon.dirname '/' oil.enveloping.name '/' oil.concat_subs.name '/' oil.concat_subs.results.concat_file_nonorm], oil.enveloping.gridstep);
else
    concat_data = nii_quickread([oil.source_recon.dirname '/' oil.enveloping.name '/' oil.concat_subs.name '/' oil.concat_subs.results.concat_file], oil.enveloping.gridstep);
end

Nsamples = size(concat_data,2);
Nics = size(spat_bas,2);
oil.ica.num_ics = Nics;
% Multiple Regression
x=demean(spat_bas,1);
pinvxtx=pinv(x'*x);
pinvx=pinv(x);
tc_est=zeros(Nics,Nsamples);
con=[]; for i=1:Nics; con{i}=(1:Nics==i)'; end  % set-up design matrix (identity).
for t=1:Nsamples;
    y=demean(concat_data(:,t));
    [tc_est(:,t)]=glm_fast_for_meg(y,x,pinvxtx,pinvx,con,0);
end

oil.ica_first_level.results.pseudo_tcs = normalise(tc_est,2);
oil.ica_first_level.results.mixing_matrix = spat_bas; 

%% Second Stage

for subnum=1:Nsubs;
    
    % Subject Specific Setup
    sing_sub_dat=nii_quickread([oil.source_recon.dirname '/' oil.enveloping.name '/' oil.enveloping.results.source_space_envelopes_NoWeightsNorm_results_fnames{oil.concat_subs.sessions_to_do(subnum)}],oil.enveloping.gridstep);
    sing_sub_tc=normalise(oil.ica_first_level.results.pseudo_tcs(:,subj_ind(subnum):subj_ind(subnum+1)-1)');  % extract subject specific part of the ICA time courses and normalise to unit std and mean.

    % Multiple Regression
    x=sing_sub_tc;
    pinvxtx=pinv(x'*x);
    pinvx=pinv(x);
    copeout=zeros(Nics,Nvoxels);varcopeout=zeros(Nics,Nvoxels);
    for v=1:Nvoxels
        y=demean(sing_sub_dat(v,:)');
        [copeout(:,v), varcopeout(:,v)]=glm_fast_for_meg(y,x,pinvxtx,pinvx,con,0);
    end
    
    % Save COPES and tstats to NII files
    cope_files{subnum} = [save_dir '/ica_copes_' oil.enveloping.results.source_space_envelopes_results_fnames{oil.concat_subs.sessions_to_do(subnum)}];
    nii_quicksave(copeout',cope_files{subnum},oil.enveloping.gridstep);
    tstat_files{subnum} = [save_dir '/ica_tstats_' oil.enveloping.results.source_space_envelopes_results_fnames{oil.concat_subs.sessions_to_do(subnum)}];
    nii_quicksave(copeout'./sqrt(varcopeout'),tstat_files{subnum},oil.enveloping.gridstep);
end

oil.ica_first_level.results.cope_files = cope_files;
oil.ica_first_level.results.tstat_files = tstat_files;
